# packages
library(tidyverse)
library(lme4)
library(psych)
library(factoextra)
library(FactoMineR)
library(corrr)
library(dotwhisker)
library(gridExtra)
library(ggpubr)
library(sjPlot)
library(sjmisc)
library(GGally)
library(mapproj)
library(rgdal)
library(plyr)
library(viridis)
library(viridisLite)
library(GPArotation)
library(ltm)
library(broom)
library(stargazer)
library(ordinal)
library(texreg)

setwd('C:/Users/test/Desktop/Lucas/Ciência Política UFPE - Mestrado/Publicações e eventos/2019/lapop/paper/scripts and data')

##### LOADING DATASETS #####
bd <- read.csv('./analysis data/ab_04-16.csv')
bd <- read.csv('../analysis data/data_reg_irt3.csv')

attach(bd)

##### DESCRIPTIVE PLOTS #####
##### MAPS #####
names(bd)

# very dissatisfied with democracy
table(very_dis)

# calculating proportion by country
map_data <- ddply(.data = bd %>% drop_na(stf_dem_num), .var = c("pais", 'wave'), .fun = function(x) {
        data.frame(n = nrow(x),
                   diss = nrow(subset(x, stf_dem_num %in% c(3, 4))),
                   diss_prop = nrow(subset(x, stf_dem_num %in% c(3, 4))) / nrow(x),
                   very_dis = nrow(subset(x, stf_dem_num == 4)),
                   very_diss_prop = nrow(subset(x, stf_dem_num == 4)) / nrow(x),
                   cor_very_wide = nrow(subset(x, bd$cor_very_wide == 1)),
                   cor_very_wide_prop = nrow(subset(x, bd$cor_very_wide == 1)) / nrow(x)
        )
}
)

bd$cor_very_wide <- ifelse(cor_p_num == 1, 1, 0)
describeBy(bd$cor_very_wide, pais)


# importing shapefile
shapefile_w <- readOGR("./original data/shapefiles/TM_WORLD_BORDERS-0.3.shp")

shapefile_w@data

# converting shapefile to dataframe 
shapefile_df <- fortify(shapefile_w)

shapefile_data <- fortify(shapefile_w@data)

# adding variable for the merge
shapefile_data$id <- row.names(shapefile_data) 

# merging geographical data
shapefile_df <- full_join(shapefile_df, shapefile_data, by="id")

# filtering countries
shapefile_df <- shapefile_df %>% filter(ISO3 %in% c("ARG", "BOL", 'BLZ', "BRA", "CHL", "COL", "CRI", "CUB", "ECU", "SLV", "GTM", "HTI", "HND", "MEX", "NIC",
                                                     "PAN", "PRY", "PER", "DOM", "URY", "VEN"))

# adding variable to the political data for merge
map_data$NAME <- map_data$pais

# merging geographic and political data
shapefile_df <- left_join(shapefile_df, map_data, by = 'NAME')

# plotting map
# dissatisfied with democracy
map1 <-  ggplot() + geom_polygon(data = shapefile_df,
                                 aes(x = long, y = lat, group = group, 
                                     fill = diss_prop * 100),
                                 colour = "lightgrey", size = .01) +
        theme_void(base_size = 11) +
        scale_fill_viridis(option = 'magma',
                           name = 'Dissatisfied (%)', direction = -1,
                           limits = c(10, 70),
                           breaks = seq(10, 70, 10),
                           guide = guide_colorbar(
                                   direction = "horizontal",
                                   barheight = unit(2, units = "mm"),
                                   barwidth = unit(50, units = "mm"),
                                   draw.ulim = F,
                                   title.position = 'top',
                                   title.hjust = 0.5,
                                   label.hjust = 0.5
                           )) +
        coord_map() +
        theme(legend.position = "bottom",
              plot.title = element_text(size = 12, face = 'bold', hjust = .5),
              plot.caption = element_text(hjust = .5))

# very dissatisfied with democracy
map2 <-  ggplot() + geom_polygon(data = shapefile_df,
                                 aes(x = long, y = lat, group = group, 
                                     fill = very_diss_prop * 100),
                                 colour = "lightgrey", size = .01) +
        theme_void(base_size = 11) +
        scale_fill_viridis(option = 'magma',
                           name = 'Very dissatisfied (%)', direction = -1,
                           guide = guide_colorbar(
                                   direction = "horizontal",
                                   barheight = unit(2, units = "mm"),
                                   barwidth = unit(50, units = "mm"),
                                   draw.ulim = F,
                                   title.position = 'top',
                                   title.hjust = 0.5,
                                   label.hjust = 0.5
                           )) +
        coord_map() +
        theme(legend.position = "bottom",
              plot.title = element_text(size = 12, face = 'bold', hjust = .5),
              plot.caption = element_text(hjust = .5))

# corruption perception in Latin America
map3 <-  ggplot() + geom_polygon(data = shapefile_df,
                                 aes(x = long, y = lat, group = group, 
                                     fill = cor_very_wide_prop),
                                 colour = "lightgrey", size = .01) +
        theme_void(base_size = 11) +
        scale_fill_viridis(option = 'magma',
                           name = 'Very widespread (%)', direction = -1,
                           guide = guide_colorbar(
                                   direction = "horizontal",
                                   barheight = unit(2, units = "mm"),
                                   barwidth = unit(50, units = "mm"),
                                   draw.ulim = F,
                                   title.position = 'top',
                                   title.hjust = 0.5,
                                   label.hjust = 0.5
                           )) +
        coord_map() +
        theme(legend.position = "bottom",
              plot.title = element_text(size = 12, face = 'bold', hjust = .5),
              plot.caption = element_text(hjust = .5))

ggsave('./command files/plots/diss_map.png', plot = map1, device = 'png',
       width = 3, height = 4.83)
ggsave('./command files/plots/ery_diss_map.png', plot = map2, device = 'png',
       width = 3, height = 4.83)
ggsave('./command files/plots/corr_map.png', plot = map3, device = 'png',
       width = 4, height = 6.44)

##### LINE GRAPHS #####

mean(map_data$diss_prop)
mean(map_data$very_diss_prop)

# very dissatisfied with democracy - average all countries
l1 <- ggplot(data = map_data, aes(x = wave, y = very_diss_prop*100)) +
        stat_summary(fun.y = mean, geom = 'line', size = 1.5) +
        theme_sjplot2() +
        labs(x = '', y = 'Very dissatisfied (%)') +
        scale_x_continuous(breaks = seq(2004, 2016, 2))

describeBy(map_data$very_diss_prop, map_data$wave)

# very dissatisfied with democracy by country
l2 <- ggplot(data = map_data %>% dplyr::filter(pais %in% c('Paraguay', 'Brazil', 'Mexico', 'Ecuador', 'Venezuela'))) +
        geom_line(aes(x = wave, y = very_diss_prop*100, colour = pais)) +
        theme_sjplot2() +
        scale_x_continuous(breaks = seq(2004, 2016, 2)) +
        scale_color_discrete(name = 'Country') +
        labs(x = '', y = 'Very dissatisfied (%)') +
        theme(legend.position = 'bottom')

l3 <- l2 +
        theme(legend.position = 'bottom') +
        geom_line(aes(x = wave, y = very_diss_prop*100, colour = pais), size = 1)

# exporting plots
#dir.create('./command files/plots')

ggsave('./command files/plots/very_diss_ctr.png', plot = l2, device = 'png', 
       width = 6, height = 4)
ggsave('./command files/plots/very_diss_mean.png', plot = l1, device = 'png',
       width = 6, height = 4)
ggsave('./command files/plots/very_diss_ctr2.png', plot = l3, device = 'png',
       width = 6, height = 4)


##### BAR PLOTS ######
table(coup_corr)
coup <- ddply(.data = bd, .var = c("pais", 'wave'), .fun = function(x) {
        data.frame(n = nrow(x),
                   coup_sup = nrow(subset(x, coup_sup == 1)),
                   coup_sup_prop = nrow(subset(x, coup_sup == 1)) / nrow(x),
                   coup_crime = nrow(subset(x, coup_crime == 1)),
                   coup_crime_prop = nrow(subset(x, coup_crime == 1)) / nrow(x),
                   coup_corr = nrow(subset(x, coup_corr == 1)),
                   coup_corr_prop = nrow(subset(x, coup_corr == 1)) / nrow(x)
        )
}
)

coup[coup == 0] <- NA

# renaming 'belice'
#coup$pais[coup$pais == 'Belice'] <- 'Belize'

# aggregating dataset by country
coup2 <- aggregate(coup, by = list(coup$pais), FUN = 'mean', na.rm = T)

# renaming variable
coup2$pais <- coup2$Group.1 %>% as.character()
coup2$Group.1 <- NULL

# proportion of supporters of military coup in case of generalized corruption
b1 <- ggplot(data = coup2 %>% drop_na(coup_corr_prop), 
             aes(x = reorder(pais, -coup_corr_prop), 
                              y = coup_corr_prop*100)) +
        geom_bar(stat = 'identity') +
        geom_hline(yintercept = base::mean(coup2$coup_corr_prop*100, na.rm = T), 
                   linetype = 2, color = 'red') +
        theme_sjplot2() +
        theme(axis.text.x = element_text(angle = 90)) +
        labs(x = '', y = '(%)') +
        scale_y_continuous(limits = c(0, 70), breaks = seq(0, 70, 10))

summary(coup2$coup_corr_prop)

# exporting
ggsave('./command files/plots/bar_coup_sup.png', plot = b1, device = 'png', 
       width = 5, height = 3.09)

# most important problem - corruption
summary(a4)

# creating dataset with proportions
coup_mip <- ddply(.data = bd %>% drop_na(a4), .var = c("pais"), .fun = function(x) {
        data.frame(n = nrow(x),
                   corruption_mip = nrow(subset(x, a4 %in% c('Corruption', 'Corrupción'))),
                   corruption_mip_prop = nrow(subset(x, a4 %in% c('Corruption', 'Corrupción'))) / nrow(x)
        )
}
)

# bar plot
b2 <- ggplot(data = coup_mip %>% filter(!pais %in% c('United States', 'CAnada')), 
             aes(x = reorder(pais, -corruption_mip_prop), 
                                  y = corruption_mip_prop*100)) +
        geom_bar(stat = 'identity') +
        geom_hline(yintercept = mean(coup_mip$corruption_mip_prop*100), linetype = 2,
                   colour = 'red') +
        labs(x = '', y = '(%)') +
        theme_sjplot2() +
        theme(axis.text.x = element_text(angle = 90))

ggsave('./command files/plots/bar_corr_mip.png', plot = b2, device = 'png',
       width = 5, height = 3.09)

###### CORRELATIONS ######

# spearman correlation between
# cor1 <- bd %>% select(cor_p_num, dstf_dem) %>%
        ggpairs(lower = 'blank',
                diag =  list(continuous = wrap(ggally_densityDiag, adjust = 7)), 
                upper = list(continuous = wrap("cor", method = 'kendall')))

# cor <- bd %>% select(cor_p_num, dstf_dem) %>% drop_na() %>%
        cor(method = 'kendall')

table(cor_p)
table(stf_dem)

##### CORRESPONDENCE ANALYSIS ######

# contingency table
tab <- bd %>% select(cor_p, stf_dem) %>% table()

# transforming to dataframe
tab2 <- as.data.frame.matrix(tab)

# correspondence analysis
ca1 <- CA(tab2, graph = T)

# plotting
ca_biplot <- fviz_ca_biplot(ca1, repel = T) +
        ggtitle('') +
        theme_minimal()

# exporting
ggsave('./command files/plots/biplot_ca.png', plot = ca_biplot, device = 'png',
       width = 6, height = 3.70)

##### FACTOR ANALYSIS #####

# scree plot
pa <- bd %>% select(b4, b43, b6, b2, b6, b11,
                    b13, b14, b21, b21a, b21e) %>% drop_na() %>%
        fa.parallel(fm = 'minres', fa = 'fa')

# extracting factors - trust in political institutions
factors <- bd %>% select(b4, b6, b2, b11,
                         b13, b14, b21, dstf_dem) %>%
        dplyr::rename(pride_pol = b4, 
                      supp_pol = b6, 
                      respect_pol = b2, 
                      trst_elec = b11, 
                      trst_leg = b13, 
                      trst_gov = b14, 
                      trst_prt = b21) %>% 
        drop_na() %>%
        fa(nfactors = 1, rotate = 'varimax')

# naming index - index of satisfaction with democracy
colnames(factors$loadings) <- 'ISD'

# correlation diagram
fa.diagram(factors)

# saving scores
bd$ISD <- factors$scores

# scree plot - only coup support
pa2 <- bd %>% 
        select(coup_sup, coup_crime, coup_corr) %>% 
        drop_na() %>%
        fa.parallel(fm = 'minres', fa = 'fa')

# extracting factors - support coup
factors2 <- bd %>% 
        select(coup_sup, coup_crime, coup_corr) %>% 
        drop_na() %>%
        fa(nfactors = 1, rotate = 'varimax')

# renaming factor - CSI: coup support index
colnames(factors2$loadings) <- 'CSI'

fa.diagram(factors2)

# two factors
factors3 <- bd %>% select(coup_crime, coup_corr, b4, b6, b2, b11,
                          b13, b21) %>%
        dplyr::rename(pride_pol = b4, 
               supp_pol = b6, 
               respect_pol = b2, 
               trst_elec = b11, 
               trst_leg = b13, 
               trst_prt = b21) %>% drop_na() %>%
        fa(nfactors = 2, rotate = 'varimax')

# renaming factors
colnames(factors3$loadings) <- c('ISD', 'CSI')

# correlation plot
fa.diagram(factors3, main = '', digits = 2, e.size = .075, rsize = 1,
           marg = c(.2, 1, .2, .2))


# saving scores in the dataset
bd2 <- bd %>% drop_na(coup_crime, coup_corr, b4, b6, b2, b11,
                      b13, b21)

bd2$CSI <- factors3$scores

# saving dataset for regression analysis
write.csv(bd2, './analysis data/data_reg.csv')

###### ITEM RESPONSE THEORY #####

# Fit graded response model for dissatisfaction with democracy
irt1 <- grm(bd[, c('b11', 'b13', 'b21')])

summary(irt1)

# Test goodness of fit (following standard cutoff of 3.5 for residuals)
margins(irt1, rule = 3.5)

# Plot item characteristic curves

plot(irt1, xlab="Theta", cex.main = .9, main = '')

# saving factors
irtscores <- factor.scores(irt1)
print(irtscores)

# getting factor values for each survey observation
ThetaTable <- factor.scores(irt1, 
                            resp.patterns = data.frame(lapply(bd[, c('b11', 'b13', 'b21')], as.numeric)))[["score.dat"]]

# binding with dataframe
bd <- cbind(bd, ThetaTable[, 7:8])

# Fit latent trait model for regime support
irt2 <- grm(bd[, c('b4', 'b6', 'b2', 'ing4')])

summary(irt2)

# Test goodness of fit (following standard cutoff of 3.5 for residuals)
margins(irt2, rule = 3.5)

# Plot item characteristic curves
plot(irt2, xlab="Theta", cex.main = .9, main = '')

# saving factors
irtscores <- factor.scores(irt2)
plot(irtscores)

# getting factor values for each survey observation
ThetaTable <- factor.scores(irt2, 
                            resp.patterns = data.frame(lapply(bd[, c('b6', 'b2', 
                                                                     'ing4')], 
                                                              as.numeric)))[["score.dat"]]

colnames(ThetaTable) <- c("b6", "b2", "ing4", "Obs", "Exp", "z3", "se.z3")

bd <- cbind(bd, ThetaTable[, 6:7])

# fit rasch model
irt3 <- rasch(bd[, c("coup_sup", 'coup_crime', 'coup_corr')])
summary(irt3)

# Test goodness of fit (following standard cutoff of 3.5 for residuals)
margins(irt3, rule = 3.5)


# Plot item characteristic curves of both models
plot(irt3, xlab="Theta")
text(x = 1.5, y = .4, labels = 'coup_sup')

# getting scores for each combination of responses
irt3scores <- factor.scores(irt3)

# extracting factors from the rasch model (the ltm model does not compute because
# of computational capacity)
ThetaTable3 <- factor.scores(irt3, 
                            resp.patterns = data.frame(lapply(bd[, c("coup_sup", 
                                                                     'coup_crime', 
                                                                     'coup_corr')], 
                                                              as.numeric)))[["score.dat"]]

# binding with main dataset
colnames(ThetaTable3) <- c('coup_sup', 'coup_crim', 'coup_corr', 'Obs', 'Exp', 
                           'z2_coup_sup', 'se.z2')

bd <- cbind(bd, ThetaTable3[, 6:7])


# exporting tables
stargazer(irt2, type = 'latex', out = '../../command files/irt_regsup.tex')

# writing dataset with IRT
write.csv(bd, '../../analysis data/data_reg_irt3.csv')

##### REGRESSIONS #####
# function to get standardized coefficients
stdCoef.merMod <- function(object) {
        sdy <- sd(getME(object,"y"))
        sdx <- apply(getME(object,"X"), 2, sd)
        sc <- fixef(object)*sdx/sdy
        se.fixef <- coef(summary(object))[,"Std. Error"]
        se <- se.fixef*sdx/sdy
        return(data.frame(stdcoef=sc, stdse=se))
}

bd$cor_p <- relevel(bd$cor_p, ref = 'Not widespread at all')
bd$soct1 <- relevel(bd$soct1, ref = 'Very Good')

bd$stf_dem <- factor(bd$stf_dem, levels = rev(levels(bd$stf_dem)))

# hierarchical ordered logit
lm1 <- clmm(data = bd, stf_dem ~ cor_p + gender + q2 + l1 + soct1 + (1 | pais) + 
                    (1 | wave))
summary(lm1)
# hierarchical logit
logit2 <- glmer(data = bd, very_dis ~ cor_p + gender + q2 + l1 + soct1 + 
                       (1 | pais) + (1 | wave), family = 'binomial')
summary(logit2)
stdCoef.merMod(logit2)

# hierarchical OLS with regime support
ols1 <- lmer(data = bd, z3 ~ cor_p + gender + q2 + l1 + soct1 + 
                     (1 | pais) + (1 | wave))
summary(ols1)
summary(bd$z3)
stdCoef.merMod(ols1)

# hierarchical OLS for authoritarianism
ols2 <- lmer(data = bd, z2_coup_sup ~ cor_p + gender + q2 + l1 + soct1 + 
                     (1 | pais) + (1 | wave))
summary(ols2)
stdCoef.merMod(ols2)

# logit with support for coup in case of widespread corruption
logit3 <- glmer(data = bd, coup_corr ~ cor_p + gender + q2 + l1 + soct1 + 
                        (1 | pais) + (1 | wave), family = 'binomial')
summary(logit3)

logit3_2 <- glmer(data = bd, coup_corr ~ cor_p_num + gender + q2 + l1 + soct1 + 
                          (1 | pais) + (1 | wave), family = 'binomial')
summary(logit3_2)

# logit with support for military coup
logit4 <- glmer(data = bd, coup_sup ~ cor_p + gender + q2 + l1 + soct1 + 
                        (1 | pais) + (1 | wave), family = 'binomial')
summary(logit4)

coef <- dwplot(list(lm1, logit2, logit3, logit4), show_intercept = F,
               vline = geom_vline(xintercept = 0), linetype = 2, colour = 'grey 60')

##### COEFFICIENT PLOTS #####
# converting model to data frames
logit1_tidy <- tidy(lm1)
logit2_tidy <- tidy(logit2)
logit3_tidy <- tidy(logit3)
logit4_tidy <- tidy(logit4)
ols1_tidy <- tidy(ols1)
ols2_tidy <- tidy(ols2)

# modifying model names
logit1_tidy <- logit1_tidy %>% mutate(model = '(1) ord logit')
logit2_tidy <- logit2_tidy %>% mutate(model = '(2) logit very diss')
logit3_tidy <- logit3_tidy %>% mutate(model = '(3) logit coup_corr')
logit4_tidy <- logit4_tidy %>% mutate(model = '(4) logit coup_sup')
ols1_tidy <- ols1_tidy %>% mutate(model = '(5) OLS regime sup')
ols2_tidy <- ols2_tidy %>% mutate(model = '(6) OLS authorit')

# adding confidence intervals
logit1_tidy$int <- confint.merMod(lm1)
logit2_tidy$int <- confint.merMod(logit2)
logit3_tidy$int <- confint.merMod(logit3)
logit4_tidy$int <- confint.merMod(logit4)

# binding models
all_models <- bind_rows(logit1_tidy, logit2_tidy, logit3_tidy, logit4_tidy)
all_models <- all_models[-c(1:3),]

# dataset with predicted probabilities for logistic regression
all_models2 <- bind_rows(logit1_tidy, logit2_tidy, logit3_tidy, logit4_tidy)
all_models2$estimate <- exp(all_models2$estimate) / (1 + exp(all_models2$estimate))

all_models2 <- all_models2[-c(1:3),]

# dataset with OLS regressions
all_models3 <- bind_rows(ols1_tidy, ols2_tidy)

# re-labeling predictors
all_models <- all_models %>% 
        relabel_predictors(c('cor_pVery widespread' = 'Corrpt wide.',
                             'cor_pSomewhat widespread' = 'Corrpt some. wide.',
                             'cor_pNot very widespread' = 'Corrpt not very wide',
                             'gender' = 'Gender',
                             'q2' = 'Age',
                             'l1' = 'Ideology',
                             'soct1Good' = 'Ecnmy good',
                             'soct1Neither good nor bad (Fair)' = 'Ecnmy fair',
                             'soct1Bad' = 'Ecnmy bad',
                             'soct1Very bad (terrible)' = 'Ecnmy terrible',
                             'sd_(Intercept).pais' = 'S.D. country',
                             'sd_(Intercept).wave' = 'S.D. Year'))

all_models2 <- all_models2 %>% 
        relabel_predictors(c('cor_pVery widespread' = 'Corrpt wide.',
                             'cor_pSomewhat widespread' = 'Corrpt some. wide.',
                             'cor_pNot very widespread' = 'Corrpt not very wide',
                             'gender' = 'Gender',
                             'q2' = 'Age',
                             'l1' = 'Ideology',
                             'soct1Good' = 'Ecnmy good',
                             'soct1Neither good nor bad (Fair)' = 'Ecnmy fair',
                             'soct1Bad' = 'Ecnmy bad',
                             'soct1Very bad (terrible)' = 'Ecnmy terrible',
                             'sd_(Intercept).pais' = 'S.D. country',
                             'sd_(Intercept).wave' = 'S.D. Year'))

all_models3 <- all_models3 %>% 
        relabel_predictors(c('cor_pVery widespread' = 'Corrpt wide.',
                             'cor_pSomewhat widespread' = 'Corrpt some. wide.',
                             'cor_pNot very widespread' = 'Corrpt not very wide',
                             'gender' = 'Gender',
                             'q2' = 'Age',
                             'l1' = 'Ideology',
                             'soct1Good' = 'Ecnmy good',
                             'soct1Neither good nor bad (Fair)' = 'Ecnmy fair',
                             'soct1Bad' = 'Ecnmy bad',
                             'soct1Very bad (terrible)' = 'Ecnmy terrible',
                             'sd_(Intercept).pais' = 'S.D. country',
                             'sd_(Intercept).wave' = 'S.D. Year')) 
table(all_models$term)

# coefficient plot 1 - log odds
coef1 <- dwplot(all_models, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = 'grey 60', linetype = 2)) +
        theme_sjplot2() +
        theme(legend.position = 'right') +
        scale_color_discrete(name = 'Models')

# coefficient plot 2 - predicted probabilities
coef2 <- dwplot(all_models2, show_intercept = F, 
                vline = geom_vline(xintercept = .5, colour = 'grey 60', linetype = 2),) +
        theme_sjplot2() +
        theme(legend.position = 'right') +
        scale_x_continuous(breaks = seq(0, 1, .25)) +
        scale_color_discrete(name = 'Models')

# coefficient plot 3 - OLS regressions        
coef3 <- dwplot(all_models3, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = 'grey 60', linetype = 2)) +
        theme_sjplot2() +
        theme(legend.position = 'right') +
        scale_color_discrete(name = 'Models')                                                    

# plotting predicted probabilities
probs1 <- plot_model(lm1, type = 'pred', terms = 'cor_p') +
        theme_sjplot2() +
        labs(x = 'Corruption perception', y = 'Satisfaction with democracy', title = '')

ggsave('./command files/plots/coef_logodds.png', plot = coef1, device = 'png', width = 6, height = 5)
ggsave('./command files/plots/coef_probs.png', plot = coef2, device = 'png', width = 6, height = 5)
ggsave('./command files/plots/coef_ols.png', plot = coef3, device = 'png', width = 6, height = 5)
ggsave('./command files/plots/probs.png', plot = probs1, device = 'png', width = 7, height = 5)



probs2 <- plot_model(logit2, type = 'pred', terms = 'cor_p_num') +
        theme_sjplot2() +
        labs(x = 'Corruption perception', y = 'Very dissatisfied (Probability)', title = '')

ggsave('./command files/plots/pred_probs.png', plot = probs2, device = 'png', width = 6, height = 4)


probs3 <- plot_model(logit3, type = 'pred', terms = 'cor_p')

probs4 <- plot_model(logit4, type = 'pred', terms = 'cor_p')

probs5 <- plot_model(ols1, type = 'pred', terms = 'cor_p')

###### MODEL TABLES #####
stargazer(lm1, logit2, logit3, logit4, title = 'Logit coefficients (log odds)',
          out = 'table_logodds.tex', model.names = T)

# table with logistic regressions
texreg(list(lm1, logit2, logit3, logit4), file = 'table_logodds.tex', bold = 0.05, 
       omit.coef = '(Very Satisfied|Very Dissatisfied) | (Very Dissatisfied|Satisfied) | 
       (Satisfied|Dissatisfied) | (Intercept)', 
       caption.above = 'Logit coefficients (log odds)', 
       groups = list('Corruption' = 1:3, 'Controls' = 4:10))

# table with OLS regressions
texreg(list(ols1, ols2), file = 'table_ols.tex', bold = 0.05, 
       caption = 'OLS coefficients', groups = list('Corruption' = 1:3,
                                                   'Controls' = 4:10))
